function [Stat,Hist]=MarkovChain_Moment(Grid,TrMat,DiffGaps)

%% Ergodic Distribution
ErgoDist    =   MarkovChain_InvDist(TrMat,[]);

%% Moment of Difference
Grid_Diff   =   bsxfun(@minus,Grid,Grid');
[qq,~,qqIdx]=   uniquetol(Grid_Diff(:),1e-6);
N_qq        =   length(qq);
N_DiffOrder =   length(DiffGaps);
Hist        =   cell(N_DiffOrder,1);
Stat        =   zeros(N_DiffOrder,4);
for ii=1:length(DiffGaps)
    if DiffGaps(ii)==0
        QQ      =   Grid;
        QW      =   ErgoDist;
    else
        QQ      =   qq;
        TempQW  =   bsxfun(@times,TrMat^DiffGaps(ii),ErgoDist');
        TempQW  =   TempQW(:);
        QW      =   zeros(N_qq,1);
        for jj=1:N_qq
            QW(jj)  =   sum(TempQW(qqIdx==jj));
        end
    end
    [QQ,TempInd]=   sort(QQ);
    QW          =   QW(TempInd);
    [CentMom,Mean,~] ...
                =   DistApp_QW_QW2Moment(QW,QQ,[1;2;3;4]);
    Std         =   sqrt(CentMom(2));
    Skew        =   CentMom(3)/Std^3;
    Kurt        =   CentMom(4)/Std^4;
    Stat(ii,:)  =   [Mean,Std^2,Skew,Kurt];
    
    Hist{ii}    =   struct('QQ',QQ,'QW',QW);
end